Discrete state space method and modal extension method based impact sound synthesis model
Cite this Article
Tian Xu-Hua, Chen Ke-An, Zhang Yan-Ni, Li Han, Xu Jian. Discrete state space method and modal extension method based impact sound synthesis model. Chinese Physics B, 2018, 27(11): 114302 
Permissions
Discrete state space method and modal extension method based impact sound synthesis model
† Corresponding author. E-mail:
Project supported by the National Natural Science Foundation of China (Grant Nos. 11574249 and 11874303) and the Natural Science Basic Research Plan in Shaanxi Province of China (Grant No. 2018JQ1001).
Abstract
The efficient and accurate synthesis of physical parameter-controllable impact sounds is essential for sound source identification. In this study, an impact sound synthesis model of a cylinder is proposed based on discrete state space (DSS) method and modal extension method (MEM). This model is comprised of the whole three processes of the physical interaction, i.e., the Hertz contact process, the transient structural response process, and the sound radiation process. Firstly, the modal expanded DSS equations of the contact system are constructed and the transient structural response of the cylinder is obtained. Then the impact sound of the cylinder is acquired using improved discrete Raleigh integral. Finally, the proposed model is verified by comparing with existing models. The results show that the proposed impact sound synthesis model is more accurate and efficient than the existing methods and easy to be extended to the impact sound synthesis of other structures.
1. Introduction
Radiated sounds in our daily life usually vary with the physical parameters of vibrating structures. These physical parameters include size, shape, and material, as well as stimuli and radiation condition.[1] Investigating the relationships between physical parameters of structures and sound features is meaningful and attractive in the field of sound source identification.[2] In order to grasp the most essential features of the vibratory and acoustical problem, sound continua should be created since they can simulate the continuous changes of the physical parameters.[3] Recordings or real sounds are not efficient to create continua since actual physical parameters are in inhomogeneous and disperse distribution.[4,5] By contrast, synthetic sounds, which can precisely control different physical parameters and provide large sets of sound samples, seem more suitable. In previous literatures, the synthesized impact sounds are mainly used in sound source identification studies since they are brodband and carrying most of the information related to physical parameters.[6,7] Thus, the synthesis of impact sound continua becomes the premised issue in the study of sound source identification.
A complete impact sound synthesis model should comprise three processes, the contact process of the hammer and the structure, the vibration response process of the structure, and the sound radiation process. The contact process was substituted by pre-defined contact forces in wave forms of a half sine or two half-Gaussians which rise fast and decay slow in some studies.[5,8] Meanwhile, the other studies used the Hertz contact model, which is closer to the reality compared with pre-defined forces, to stimulate the whole processes of the contact.[9,10] The vibration response and radiation process was usually modeled with basic theoretical derivation, modal expansion method (MEM), finite difference method (FDM), finite element method (FEM), boundary element method (BEM), and combination of the above methods.[9–13] Lambourg (2001) synthesized the vibration response and impact sound through the finite-difference time-domain (FDTD) method and discrete Raleigh integral, but his method had a drawback of low computational efficiency.[9] Lutfi (2001) and McAdams (2010) used simplified models based on Lambourg’s method to synthesize the impact sound of bar and plate for subjective judgments, but their models would not be suitable for the accurate-sound-feature extraction since they ignored the information which is not sensitive by the human ears.[10,11] Zhang (2014) proposed FDTD-MEM method and FEM-MEM-BEM method (implemented in LMS Virtual.Lab software) to synthesize impact sound more efficiently than FDTD method under the same precision.[12,13] However, Zhang’s methods still suffered from low computational efficiency caused by high order FDTD and large number of finite elements. In addition, Avanzini (2001) put forward a low parameter sound synthesis model, from which the system was constructed by DSS equations discretized through bilinear transform method.[14] The DSS equations are first order with lower computing cost than high order FDTD equations. Nevertheless, Avanzini’s model neglects both stimulus-position information and radiation process, and the model is not accurate due to the bilinear transform method adopted to discrete the state space equations. Although aforementioned methods have advantages respectively in accuracy and efficiency, an efficient and accurate method for sound synthesis still needs to be developed for sound feature identification tasks.
In order to have a simultaneous improvement on the accuracy and the computational efficiency compared to aforementioned methods and ensure the information integrity of the synthesized sounds, an impact sound synthesis model is proposed in this study. The model comprises the Hertz contact process, the vibration response process, and the sound radiation process. And it is firstly built with modal expanded DSS equations. Then the vibration response and radiated sound of the simply-supported cylinder, which is widely used in airplane and ship manufactures, are calculated using two approximate methods (Euler method and trapezoid method). The results show that the proposed sound synthesis model is superior to the existing methods in both accuracy and efficiency. The proposed model can provide a large number of sound samples for the subsequent sound source identification study and it can be extended to other structures with just a few modifications of the modal parameters.
2. Impact sound synthesis model
2.1. Overview of the model
The impact sound synthesis model used in this paper is described briefly in Fig.
 | Fig. 1. The sound synthesis model. |
2.2. Basic equations
Basic equations of different parts of the sound synthesis model are displayed in the following subsections. They will be discretized and solved in a form of regressive iteration equations in Section
2.2.1. Cylinder motion
where the superscript (r) represents the resonant cylinder,
,
,
, and
are respectively the modal amplitude, the modal damping, the modal angular frequency, and the modal forces of the l-th mode. By defining varphi l as the l-th modal shape of the structure, the terms in Eq. (1 ) can be described as follows:[15]
where km, n, θ, ρ, and h stand for the axial wave number, the circumferential order, the circumferential angle, the mass density, and the thickness of the cylinder, respectively. ζl denotes the l-th (corresponding to the (m,n) order mode of the cylinder) modal damping ratio. The total displacement w(xp,t) of point xp can be calculated according to the modal analysis theory as
where φl(xp) is the positional coefficient of the l-th modal shape at point xp.
In order to decrease the dimension of iterative computation, the equations of the contact system are usually modal expanded using MEM. The modal expanded second order differential equation of the cylinder can be expressed as follows:
 |
 |
 |
2.2.2. Hammer motion
where the superscript h represents the hammer, and f(t) is the contact force.
Assuming the contact process is on the horizontal direction, the equation of hammer motion can be described as
 |
2.2.3. Contact force
where k is the contact stiffness, λ and α denote the damping coefficient and the index, respectively, and they are set to 0.6k and 3.8 following Avanzini’s configurations. Δx(t) and Δv(t) represent the relative displacement and the relative velocity between the hammer and the cylinder respectively
Contact force (or impact force) is one of the crucial clues associated with the physical parameters of the hammer, such as the hardness.[16] As mentioned in Section
 |
3. Structural response based on DSS and MEM
3.1. Theory
where
Instead of discretizing the equations directly as Avanzini did,[14] we solve Eq. (8 ) and acquire an accurate solution before the discretization operation. The accurate solution of Eq. (8 ) is shown as follows:[19]
where the first term on the right hand side is the solution of the free vibration, and the second term is the forced vibration solution (also called Duhamel integral), and τ ∈ [t0, t] is the period considered. The discrete form of Eq. (9 ) can be found in Ref. [19]
where
is the average value in τ ∈ [nT − T,nT], Fs is the sampling rate, and Θ equals to
. The upper limit of i can be set to 5–10 since
is negligible when i exceeds 5. To improve the iterative efficiency, equation (10 ) is divided into two iterative formulas
where X (n) and V (n) are respectively the displacement amplitude vector and the velocity amplitude vector of all the modes considered. eA ′T and Θ ′B ′/Fs are pre-calculated before the iteration process since they are time independent. The quantities A ′, B ′, and Θ′ are expressed as follows:
where I ′ and (A ′)i should be redefined as
where “ ° ” is the Hadamard product, i.e., the product of elements at the same position of the matrixes.
where
can be approximated by the Euler method (
or the trapezoid method (
. In general, the trapezoid method has a better accuracy in approaching the velocity gradient than the Euler method. However, the trapezoid method is more complex since u (n) and x r)(n) in Eqs. (11 ) and (16 ) have instantaneous mutual dependence. Fortunately, K method proposed by Borin[20] could be introduced to calculate u (n) in each time step with a few steps of Newton–Raphson iterations. Detail information about K method used in this study can be found in Ref. [15].
where TpAl is the transfer coefficient between different points. Thus, the l-th structural response of the whole cylinder can be calculated through Eq. (18 ). The driven point must be on the nodal line of some modes, which may cause unexpected high value of 1/φpl. Therefore, TpAl is weighed with a coefficient which equals to 0 when φpl is smaller than a certain threshold (0.01, for the proposed model) and 1 otherwise.
In order to improve the accuracy and the computational efficiency of aforementioned methods, we propose an accurate DSS and MEM based method to calculate the structural response of the cylinder. The modal expanded DSS equation of Eq. (
 |
 |
 |
 |
 |
 |
 |
 |
 |
Accordingly, the DSS equations of the hammer and the contact model can be described as
 |
 |
By solving Eqs. (
 |
3.2. Results and validation
The parameters of the contact system are given in Table
 | Table 1. Physical parameters of the system. . |
3.2.1. Results
The temporal contact force and velocity response of the cylinder from different approximate methods are displayed in Figs.
 | Fig. 2. Temporal contact force. |
 | Fig. 3. (color online) Velocity response. |
3.2.2. Validation
The driven-point admittance results of different methods are displayed in Fig. 4 in different line styles. Among them, the results of the proposed methods and bilinear transform method are calculated using Eq. (19 ), meanwhile the theoretical result is calculated from Eq. (20 ). It is noticed that the driven-point position is considered in all the above methods. The results show that the results of the proposed methods well agree with the theoretical result and the bilinear transform result has obvious deviation in both amplitude and frequency.
The FEA software can hardly work out a result under the calculation precision of 5000 Hz since the required element number for the cylinder in Table
 |
 |
 | Fig. 4. (color online) Driven-point admittances of different methods. |
The significant deviation of Avanzini's bilinear transform result can be interpreted from two aspects: 1) the Duhamel integral in Eq. (
As a complementary validation, the velocity response of a smaller cylinder (L = 0.3 R=0.1) calculated using FEM+MEM (realized through the transient modal superposition structural response case in Virtual.Lab software) and proposed method is shown in Fig.
 | Fig. 5. Contact point velocity of different methods. |
3.3. Computational efficiency
The time consumed in different steps of the transient structural response process is displayed in Table
| Table 2. Time consumed in different steps. . |
It is noticed that the proposed methods, especially the trapezoid approximate method, are more efficient than the bilinear transform method used by Avanzini. This phenomenon reveals that pre-calculation of two terms in Eq. (
4. Sound pressure calculation and comparison
4.1. Sound pressure calculation
where Pj is the sound pressure at xj, V and T (a) are respectively the normal velocity vector of discrete surface nodes and the acoustic transfer vector. V and T (a) can be described as follows;
where
is the distance between the i-th structural surface node and the j-th sound field point, ΔS is the area of single surface element, U(n) is the Heaviside function, and Round(x) rounds x to the nearest integer. The time delays of different elements are considered in Eq. (22 ). Suppose that the positions of a structural node and a sound field point are respectively (r0, θi, zi) and (r, θj, zj), the distance
can be obtained by
A classical way to calculate the sound field of a cylinder is the continuous Raleigh integral method in cylindrical coordinate,[15] but it is too complex to meet the requirement of high computational efficiency. On the contrary, discrete Raleigh integral is more efficient since the computational dimension is greatly reduced by the discretization operation.[9,13] There are two discrete forms of Raleigh integrals in Eq. (
 |
 |
 |
4.2. Result and validation
The receiving point of the sound field is set to xrecieve(r, θ,z) = (2,0,0.08) and the impact sounds of different Raleigh integral forms are shown in Fig.
 | Fig. 6. Comparison of impact sounds from different Raleigh integral forms. |
It is obvious that the impact sounds from different methods well agree with each other. The proposed first order method has the same accuracy to the second order FDTD method used by Lambourg and Zhang.
For further validation of the proposed impact sound synthesis method, the impact sound of a smaller cylinder (L = 0.3 R = 0.1) calculated respectively using FEM+MEM+BEM (realized through the transient modal superposition weakly coupled response case in Virtual.Lab software) and proposed method are shown in Fig.
 | Fig. 7. Comparison of impact sounds from different impact sound synthesis methods. |
4.3. Computation efficiency
The computing time of the impact sound (does not contain the Hertz contact and velocity response processes) is displayed in Table
| Table 3. Comparison of the impact sound computation efficiency. . |
5. Conclusion
Based on a simply supported cylindrical shell structure, an impact sound synthesis model comprising the Hertz contact process, the transient structural response process, and the sound radiation process is proposed. Firstly, the system equations are constructed through the modal extended DSS equations and the transient structural responses of the cylinder are obtained. Then the impact sound is calculated using the proposed first order discrete Rayleigh integral method. Finally, the sound synthesis model is verified in two steps: 1) verifying the structural response via the point admittance and FEM+MEM method, and 2) verifying the sound radiation by comparing impact sounds from different methods.
The results show that the DSS and MEM based sound synthesis model has advantages over the existing sound synthesis models in terms of both computational efficiency and accuracy. The sound synthesis model could provide a large number of effective sound samples for sound source identification study and could be easily transferred to the impact sound synthesis of other structures.
Reference